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Abstract. This paper is concerned with the Einstein equations in axisymmetric 
vacuum spacetimes. We consider numerical evolution schemes that solve the 
constraint equations as well as elliptic gauge conditions at each time step. We 
examine two such schemes that have been proposed in the literature and show 
that some of their elliptic equations are indefinite, thus potentially admitting 
nonunique solutions and causing numerical solvers based on classical relaxation 
methods to fail. A new scheme is then presented that does not suffer from 
these problems. We use our numerical implementation to study the gravitational 
collapse of Brill waves. A highly prolate wave is shown to form a black hole rather 
than a naked singularity. 



PACS numbers: 04.25. D-, 04.20.Cv, 04.20.Ex, 04.25.dc 

1. Introduction 

Driven by the need of gravitational wave data analysis for waveform templates, 
numerical relativity has focused in recent years on the modelling of astrophysical 
sources of gravitational waves such as the inspiral and coalescence of compact objects. 
Such systems do not possess any symmetries and thus require a fully 3-1-1 dimensional 
numerical code. The advantage of assuming a spacetime symmetry, on the other 
hand, is that it allows for a dimensional reduction of the Einstein equations, which 
reduces the computational effort considerably so that greater numerical accuracy can 
be obtained. While spherical or planar symmetry yields the greatest reduction in 
computational cost, the intermediate case of axisymmetry is more interesting in that 
it permits the study of gravitational waves. 

In this article we focus on vacuum axisymmetric spacetimes and assume that 
the Killing vector is hypersurface orthogonal so that there is only one gravitational 
degree of freedom. The axisymmetric Einstein equations can be simplified considerably 
by choosing suitable gauge conditions. Here we consider a combination of maximal 
slicing and quasi-isotropic gauge [l][2l[3j. This gauge reduces the number of dependent 
variables to such an extent that only one pair of evolution equations corresponding 
to the one gravitational degree of freedom needs to be kept. All the other variables 
can be solved for using the constraint equations and gauge conditions. This fully 
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constrained approach was taken in |4j. Partially constrained schemes (e.g., [SJ [B]) 
substitute some of the constraint equations with evolution equations; this is possible 
because the Einstein equations are overdetermined. 

Such (fully or partially) constrained schemes have proven very robust in 
simulations of strong gravity phenomena. Examples include the collapse of vacuum 
axisymmetric gravitational waves, so-called Brill waves [7], in [6]. Critical phenomena 
in this system were found in [5l[8] . Critical phenomena in the collapse of massless scalar 
fields [9] and complex scalar fields with angular momentum [10] were also studied, as 
was the collapse of collisionless matter [TT| . 

Nevertheless, constrained evolution schemes have been plagued with problems. 
The authors of [4] reported that their multigrid elliptic solver failed occasionally for 
the Hamiltonian constraint equation in the strong-field regime. This problem can 
be circumvented by using instead the evolution equation for the conformal factor. 
However, it was found that this was "not sufficient to ensure convergence in certain 
Brill- wave dominated spacetimes" . Similar difficulties were encountered in [T^ [T3] . 
The purpose of this article is to determine the cause of these problems and to develop 
an improved constrained evolution scheme. 

The suspect elliptic equations belong to a class of (nonlinear) Helmholtz-like 
equations, which are discussed quite generally in section [2l We point out that if these 
Helmholtz equations are indefinite (loosely speaking, they have the "wrong sign" ) then 
their solutions, should they exist, are potentially nonunique. The same criterion is 
found to be related to the convergence of numerical solvers based on classical relaxation 
methods. 

In section [3l we review the partially constrained scheme of [6] and the fully 
constrained scheme of 14[. We show that some of the elliptic equations in these 
formulations are indefinite. This leads us to the construction of a modified fully 
constrained scheme that does not suffer from this problem. The arguments involved 
turn out to be closely related to questions of (non)uniqueness in conformal approaches 
to the initial data problem in standard 3 -t- 1 numerical relativity [111 [15] . 

A numerical implementation of the new fully constrained scheme is described in 
section H) In section [S] we apply it to a study of Brill wave gravitational collapse. 
After performing a convergence test and comparing results for a strong wave with 
"spherical" initial data, we focus on a highly prolate configuration-one of the initial 
data sets examined in [16]. By considering sufficiently prolate configurations, the 
authors were able to construct initial data without an apparent horizon but apparently 
with arbitrarily large curvature. They conjectured that such initial data would evolve 
to form a naked singularity rather than a black hole. This would constitute a violation 
of weak cosmic censorship. A numerical evolution of one of these prolate initial data 
sets was carried out in [6] . Due to a lack of resolution on their compactified spatial 
grid, the authors could not evolve the wave for a sufficiently long time. The trends in 
certain quantities suggested however that an apparent horizon would eventually form. 
Using our new constrained evolution scheme, we are able to evolve the same initial 
data for much longer and we confirm that an apparent horizon does form. 

We conclude and discuss some open questions in section [6l 
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2. General remarks on Helmholtz-like elliptic equations 

Let us first consider the Helmholtz equation 

Au + cu^f, (1) 

where u £ M", A is the fiat-space Laplacian, and c and / are smooth functions. We 
impose the boundary condition w ^ at spatial infinity. More generally, a boundary 
condition u —i- a — const can always be transformed to this case by considering the 
function u — a. It follows from standard elliptic theory (see e.g. [17]) that ([1]) has a 
unique solution if c ^ everywhere. If c > then multiple solutions may exist or 
there may not be any solution at all. For c > the elliptic equation is said to be 
indefinite. 

Next we consider the quasilinear equation 

Au + F{u) = f, (2) 

where F(u) is a smooth (not necessarily linear) function. Proving existence and 
uniqueness of solutions to this equation is nontrivial. However a necessary condition 
for uniqueness can easily be obtained. Suppose uq is a given solution and consider a 
small perturbation of it, u — uq + Su. Approximating F{u) ~ F{uo) + F'{uo)Su where 
F'{u) = dF{u)/du, we find that for u to be a solution Su must satisfy 

Adu + F'{ua)6u^0. (3) 

This is of the form ^ with c = F'{uo) and / 0. If F'(wo) < then there is only 
the trivial solution 6u = and we call the original problem ([2]) linearization stable. 
If on the other hand F'{uo) > then multiple solutions of the linearized problem ([3]) 
and hence also of the nonlinear problem ^ may exist. 

As an example relevant to the formulations of the Einstein equations discussed 
in this article, we take 

F{u) = cuP (4) 

with p G M and c a smooth function. Then ([2]) is linearization stable provided that 
pc ^ 0. We say that in this case the equation has the "right sign" . 

Because of the above considerations on the uniqueness of solutions, it is clearly 
desirable to have an equation with the "right sign" if a numerical solution is attempted. 
There is however also a more practical reason. Consider again the linear Helmholtz 
equation ([1]) in, say, n = 2 dimensions. Suppose we cover the domain with a 
uniform Cartesian grid with spacing /i, denoting the value of u at the grid point 
with indices {i,j) by Uij. A discretization of ([T]) using second-order accurate centred 
finite differences yields 

Ui+ij + Wi-ij -t- Uij+i + Uij-1 - (4 - ch'^)uij = h^fij. (5) 

We formally write this system of linear equations as 

Au = /i^f, u = {wy}. (6) 

Large systems ([B]) are commonly solved using relaxation methods, which obtain a series 
of successively improved numerical approximations. For example, a step of the Gauss- 
Seidel method consists in sweeping through the grid (typically in lexicographical or 
in red-black order), at each grid point (i, j) solving the equation for Uij and replacing 
its value, 

Uij ^ (ui+ij + Wi-i J + Uij+i + Uij-i - h?fij)/ (4 - ch^). (7) 
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The relaxation converges provided the matrix A in is strictly diagonally dominant, 
i.e., in each row of the matrix the absolute value of the diagonal term is greater 
than the sum of the absolute values of the off-diagonal terms (see e.g. [H]). In our 
example (O, the diagonal term is |4 — c/i^| and the off-diagonal terms add up to 4, 
so that the condition for convergence is c < 0. (The other possibilty, c large and 
positive, is not feasible because /i — > in the continuum limit.) In practice, if c 
is positive but sufhciently small then the relaxation will still converge but as c is 
increased convergence will begin to stall and ultimately the relaxation will diverge. 
Similar convergence criteria hold for other relaxation schemes such as the Jacobi or 
SOR methods. In particular, the multigrid method [19] is based on these relaxation 
schemes and will not converge if the underlying relaxation does not. These warnings 
do not apply to certain versions of the conjugate gradient method [Hj or other Krylov 
subspace iterations which ideally only require the matrix A to be invertible. For a 
combination of such methods with multigrid see e.g. |20| . 

3. Formulations of the axisymmetric Einstein equations 

We focus on axisymmetric vacuum spacetimes. Axisymmetry means that there is 
an everywhere spacelike Killing vector field ^ with closed orbits. Here we restrict 
ourselves to the case where the Killing vector is hypersurface orthogonal. We choose 
cylindrical polar coordinates t,z,r,(j> such that ^ = d/d(f). In the following, indices 
a,b, . . . range over t, r, z, 4>, indices i, j, . . . over r, z, 0, and indices A,B, . . . over r, z. 
The line element is written in the form 

ds^ = -a^dt^ -f V'^{e2''-5[(dr - P^'dtf + (dz - P^'dtf] + r^d^^}. (8) 

Here a and (3^ are the usual ADM lapse function and shift vector. 

We have imposed as a gauge condition that the 2-metric on the t = const, 4> = 
const hypersurfaces be conformally flat in our coordinates (quasi-isotropic gauge) : the 
spatial metric obeys ^rz — and 7rr — Izz- This condition must be preserved by the 
evolution equation for the spatial metric, 

dtli] = ~2aKij +C(3-fij, (9) 

where Kij is the extrinsic curvature of the t = const surfaces and C denotes the Lie 
derivative. We deduce that 

f]+ = f3\r+l3'^,z = 2aK%, 

/3„ = - I3\z = ~aU, (10) 

where U = - K\. 

Maximal slicing K = Ki^ = is imposed, so that the extrinsic curvature has 
three degrees of freedom, which are taken to be K^r^ U, and W = (K^r — K'^^)/r 
(this particular combination is motivated by regularity on the axis of symmetry [6 ) . 
The evolution equation for the extrinsic curvature is given by 

dtK\j = -DWja + aR}j + CpK'-j, (11) 

where D is the covariant derivative compatible with the spatial metric 7ij and Rij is 
its Ricci tensor. Preservation of the maximal slicing condition implies DiD^a = aR = 
aK^jK^i (using the Hamiltonian constraint for the second equality) or 

a,rr + a,.;, -f- {2Pr + r-'^)a^r + 2P,a,^ 

-2a^/>V''^ [\{U + \rWf + \{rWf + {K\f] = 0. (12) 
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Here and in the following we use the notation 

Pa = ^-^^,a. RA = {rS)^A, AA = a-^aA. (13) 

There are many different ways of constructing an evolution scheme for the 
axisymmetric Einstein equations in the above gauge, depending on the number of 
constraint equations being solved. We review two schemes that have been used for 
numerical simulations and show that some of their elliptic equations are indefinite as 
discussed in the section [5] Finally we propose a new scheme that does not suffer from 
this problem. 

3.1. A partially constrained scheme: Garfinkle and Duncan JBj/ 

Garfinkle and Duncan [6] choose to solve only the Hamiltonian constraint equation 
R — jK^ i — (note that K = {)), which takes the form 

+ 1^5g2rS [1 ^ 1^^)2 ^ ^[rWf + {K\f] = 0. (14) 

This equation is of the type ^ with p = 5 and c ^ (note the second square 
bracket in (fH)) is non- negative) . Hence it has the "wrong sign" and suffers from 
potential nonuniqueness of solutions as well as difficulties in solving it numerically 
using relaxation methods (section[2]). The latter is not a concern in [6] though because 
the authors use a conjugate gradient method. 

The momentum constraints Kij = are not solved but only monitored during 
the evolution. Written out explicitly they are 

2K\^^ - lU,r + IrW^r + (12P^ + AR^)K\ - [APr + 2Rr)U + {ArPr + l)W = 0, 
2K''r.r + + IrW,^ + (12P,. + ARr + 2r-^)K\ + {8P^ + 2R^)U + AP^rW = 0. 

(15) 

The extrinsic curvature variables K^'r, U and W are all evolved using their time 
evolution equations [6]. 

The slicing condition is solved in the form and this is a Helmholtz equation 
with the "right sign" (c ^ in ([1]); note the square bracket in (fT^ is non-negative). 

In order to solve for the shift vector, additional derivatives are taken of equations 
(jlOp. which combine to two decoupled second-order equations, 

P\rr+P\z.^2{aK\),,-{aU)^r, 

P\rr + P\zz = 2{aK\)^r + ("C/),,. (16) 

These are Poisson equations (c = in ([1])) and do not cause any problems. 

3.2. A fully constrained scheme: Choptuik et al 

A similar formulation was developed by Choptuik et al [4j. Their definition of the 
variables a and ip differs slightly from our S and tp, 

a = -5, VSh = ^'e^'^^ (17) 

where the subscript 'Ch' refers to [J]. This difference does not have any consequences 
on the properties of the elliptic equations that we are concerned with here and so for 
the sake of consistency we continue to use our convention (which agrees with the one 
in [Hj). As a result the equations displayed below differ from those in [J in a minor 
way. 
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In the same way as Garfinkle and Duncan, Choptuik et al also solve the 
Hamiltonian constraint, which is again indefinite. 

Unlike Garfinkle and Duncan, however, they also solve the momentum 
constraints. This is done by replacing K^r and U with first derivatives of the shift 
using the gauge conditions (|10p . The momentum constraints (I15|) now read 

p\rr + f3\zz + \(i\rz + ^arW + a{ArPr + l)W 

+ (6P, + 2Rz - Az)l3+ + [APr + 2Rr - |A,)/3_ ^ 0, 

+ {6Pr + 2Rr + r~i - Ar)(3+ - (8P, + 2Rz - \Az)Si- = 0. (18) 

The principal part of these two coupled equations is elliptic and so far there is no 
need for concern. A problem arises however when equations (|10p are substituted in 
the shcing condition (fT2|) . 

+ a,.. + (2P,. + r-^)oc.r + 2P,a,^ - V^V'^^a-^f/^- + 

^\^^(?^^ (i-rW - iVV'^^a(rM^)2 = 0. (19) 
The term containing the square bracket has the "wrong sign" , p = — 1 and c ^ in 

3.3. A new fully constrained scheme 

We observed that in both of the above schemes, the Hamiltonian constraint was 
indefinite, and in the second one, the slicing condition was, too. We now present 
a scheme in which both equations and in fact all the elliptic equations that are being 
solved are definite. 

The Hamiltonian constraint can be cured by rescaling the extrinsic curvature 
variables with a suitable power of the conformal factor, 

{K%,U,W} = ij''{K%,U,W}. (20) 

In terms of the new variables, the exponent of ip multiplying the second square bracket 
in is 5 — 2p so that the equation becomes definite for p ^ 5/2. There is a 
preferred choice: for p = 6 the terms containing derivatives of ip in the momentum 
constraints pS]) all cancel under the substitution (PU)) . The same rescaling of the 
extrinsic curvature was applied by Abrahams and Evans 5j . Their scheme is however 
not fully constrained-the extrinsic curvature variables are evolved as in ^ . 

The indefiniteness of the slicing condition was caused by the substitution (fTO]) . 
more precisely by its a dependence. The original motivation for this substitution was 
the desire to be able to solve the momentum constraints. However we can still do this 
as before if we introduce a new vector r/^ and set 

T]+=ri\r + v\z^2K%, 

ij_ = i^\r-v\. = -U. (21) 

The momentum constraints are then solved for 77"*. The price we have to pay is that 
we still need to solve the spatial gauge conditions ((T6)) . where now K^r and U are 
expressed in terms of ri±. That is, we have to solve two more elliptic equations than 
Choptuik et al . 
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Let us now write out all the elliptic equations explicitly. The momentuni 
constraints are 

h'^rr + V\zz + + 2i?,?7_ + 2R,T^+ + frW, + f W = 0, 

-\l\rz + ?7^rr + ^V^rr - 2R,V- + {^Rr + r-^)v+ + ^rWz = 0. (22) 

The Hamiltonian constraint is 

i>,rr + tp.zz + r~^1p,r + J^[{rS)^rr + {rS),zz] 

+ 3lV'"V'''5[(2?7_ - rW)^ + 3{rW)^ + 3ry^] ^ 0. (23) 
The slicing condition is 

a,rr + azz + {2Pr + r^^)a,r + 2Pza,z 

-iai/'"*^e2'-'5[(27?_ - rWf + Z{rWY + a??^,] = 0. (24) 
The spatial gauge conditions are 

I3\rr + I3\.. - aij'%\rr + v\zz " (6P. - A,)r;_ - (6P, - A,)r/+] = 0, 

P\rr + r - aV^^'^try'.rr + ?]",22 + (6P. - A,)t]- ~ (6P, - A)?y+] = 0. (25) 

We note that form a hierarchy: the equations are successively solved for 

Ty'^, ip^ ^iid (3^. After substituting the solutions of the previous equations, each 
equation in the hierarchy can be regarded as a decoupled scalar elliptic equation, or 
elliptic system in the case of ()22p . The terms in the second lines of (|23p and ([24| 
now have the "right signs" . An exception common to all the schemes discussed in this 
section is the term multiplying -0 in the first line of (|23p -in general one expects S to 
oscillate so that A(r5) can have either sign. This is the usual difficulty one faces in 
conformal formulations of the initial value equations, see section [3.41 

The variable S and its "time derivative" W are evolved. This pair of evolution 
equations corresponds to the one dynamical degree of freedom. Note that if we had 
not restricted the Killing vector to be hypersurface orthogonal then there would be 
a second dynamical degree of freedom. In linearized theory these two degrees of 
freedom can be understood as the two polarization states of a gravitational wave. 
There are additional evolution equations for ij:, 77+ = 2K^r and rj- = —U that are not 
actively enforced but that can be used in order to test the accuracy of a numerical 
implementation. All the evolution equations are given in |Appcndix A| Here we remark 
that assuming a solution to the elliptic equations is given, the principal part of the 
evolution equations is that of a wave equation, 

[a-\dt - f3^dA)Ys ~ ^-^e-^^^iS,rr + 5,,), (26) 

where ~ denotes equality to principal parts. This equation is clearly hyperbolic, a 
necessary criterion for the well posedness of the Cauchy problem. See also \21\ for a 
recent analysis of the hyperbolic part of a fully 3 + 1 dimensional constrained evolution 
scheme based on the Dirac gauge. 

3.4- Relation to the (extended) conformal thin sandwich formulation 

Our discussion of different constrained evolution schemes for the axisymmetric 
Einstein equations is closely related to conformal approaches to solving the initial 
value equations in standard 3 + 1 dimensional spacetime. Here one seeks to find a 
spatial metric 7^^ and extrinsic curvature Kij satisfying the constraint equations on the 
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initial slice, and often also a lapse a and shift satisfying suitable gauge conditions. 
This is done by setting 

lij = i^^lij , (27) 

where ip is the conformal factor, and the conformal metric "fij is assumed to be given. 
For simplicity and for analogy with the axisymmetric formulations discussed above, we 
impose the gauge condition dtjij = 0. (In the axisymmetric case, we only controlled 
the r, z components, dtjAB — 0.) Wc also assume maximal slicing K = throughout, 
and we work in vacuum. 

It is well known in the conformal approach that the extrinsic curvature Kij cannot 
be freely specified [25 ; instead it has to be conformally rescaled, 

^y-K7y =V'"'%-- (28) 

This corresponds to the proposed rescaling of the extrinsic curvature variables (pO| in 
the new axisymmetric scheme. The Hamiltonian constraint now takes the form of the 
Lichnerowicz equation, 

BAi/' - R^J + ij-^AijA'^ = 0, (29) 

where A is the covariant Laplacian and R the Ricci scalar of the conformal metric "fij . 
Note again that the last term in ([29]) has the "right sign" for linearization stability, 
cf. ([23]) . As pointed out in the previous subsection, the linear term —Rtp can have 
either sign. However, i? < does not necessarily imply that multiple solutions exist. 
There is a well-developed theory for existence and uniqueness of solutions to ([29]) , see 



In order to solve the momentum constraints, York's original conformal transverse 
traceless (CTT) method introduces a vector rj^ and sets 

i„=(Lr/)„-, (30) 
where L is the conformal Killing operator defined as 

ihrj),^ = 2V(,?7j) - |7.jVfe?7^ (31) 
The momentum constraints now read 

V.ihfjr^O, (32) 

in analogy with ([22]) . 

In the CTT approach, any gauge conditions are solved after a solution of the 
constraint equations has been found. For example, maximal slicing K — dtK — 
implies the following elliptic equation for the conformal lapse a — ip^^a, 

Aa - AijA'^a = 0, (33) 

where ([30| is substituted. This equation has the "right sign", as it has in the new 
axisymmetric formulation (|24p and in the one by Garfinkle and Duncan (jl2p . 

In contrast, the extended conformal thin sandwich (XCTS) method [24j directly 
expresses the extrinsic curvature in terms of the shift 

4- ^ Mi (34) 

instead of ([50)1 . As a result, the slicing condition ([55]) acquires the wrong sign. This 
is precisely what happens in the scheme by Choptuik et al , cf. HUj) and (|19p . 

Remarkably, a numerical study of the XCTS equations [Mj showed that this 
system does admit nonunique solutions. Two solutions were found for small 
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perturbations of Minkowski space, one of them even containing a black hole. The 
two branches meet for a certain critical amplitude of the perturbation. This 
parabolic branching was explained using Lyapunov-Schmidt theory in [15j . Because 
of the similarity with the XCTS equations, it is conceivable that the constrained 
axisymmetric formulation of Choptuik et al [4] might show a similar branching 
behaviour. This is clearly undesirable for numerical evolutions because the elliptic 
solver might jump from one solution branch to the other during the course of an 
evolution. However even before this can happen the multigrid method used in [4J will 
fail to converge, as explained in section O 

4. Numerical method 

In this section we describe a numerical implementation of the new fully constrained 
scheme presented in section [531 

The equations are discretized using second-order accurate finite differences in 
space. A collection of the finite difference operators we use can be found in 
[Appendix B[ Similarly to 6 , and unlike [4], we use a cell-centred grid to cover the 
spatial domain [0, rmax] x [0, ^max]: grid points are placed at coordinates = (i— |)Ar, 
1 ^ J ^ Nr, where Ar — r^^^^/Nr is the grid spacing and N^. is the number of grid 
points in the r direction. (Corresponding relations hold in the z direction.) Note that 
no grid points are placed at the boundaries. Ghost points are placed at — — ^Ar 
and at rj^^+i = rmax + ^Ar. The values at these ghost points are set according to 
the boundary conditions, as described in the following. Here we only refer to the 
"physical" grid boundaries at r = 0, z = 0, r = r^ax and z = ^max- In the adaptive 
mesh refinement approach discussed further below, additional finer grids are added 
that do not cover the entire spatial domain. These finer grids are also surrounded by 
ghost points. On grid boundaries that do not coincide with a "physical" boundary, 
ghost point values are interpolated from the coarser grid [25[ . 

The boundary conditions at r = follow from regularity on axis (see [26| for a 
rigorous discussion) : either a Dirichlet or a Neumann condition is imposed depending 
on whether the variable is an odd or even function of r. All the equations being 
solved (both the elliptic equations and the evolution equations) are regular on the 
axis provided that these boundary conditions are satisfied. In addition, we impose 
reflection symmetry about z = so that the variables are either odd or even functions 
of z, and this implies Dirichlet or Neumann conditions at z = 0. The r and z parities 
of all the variables are listed in [Appendix B[ 

At the outer boundaries r — rmax and z — Zmax, we impose Dirichlet conditions 
on the gauge variables, a = 1 and /3"^ — 0. For the variables {V', V^} ^ u, we impose 

M = Uoo + ^, (35) 

where R = Vr^ + and 9 — tan(r/z) are spherical polar coordinates and Uao is the 
value of u at spacelike infinity, i.e., "000 = 1 and rj^ = 0. This boundary condition 
obviously holds up to terms of 0{R)~^ for any asymptotically flat solution of the 
constraint equations. For the "dynamical" fields {S, W} 3 u, we follow 4J and impose 
a Sommerfeld condition at the outer boundary, 

{dt + dR)[Riu - uo.)] ^ 0. (36) 

This condition is only exact for the scalar wave equation in spherical symmetry. It 
is however expected to be a reasonable first approximation because S and W obey a 
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wave equation (|26p to principal parts and the elliptic variables will be close to their 
flat-space values near the outer boundary {a — ifi — 1, [3^ — -q^ = 0). Sec [Appendix B| 
for details on the discretization at the outer boundary. 

The evolution equations are integrated forward in time using the Crank-Nicholson 
method; this method is second-order accurate in time. The resulting implicit equations 
are solved by an outer Gauss-Seidel relaxation (in red-black ordering), and an inner 
Newton-Raphson method in order to solve for the vector of unknowns at each grid 
point. A typical value of the CFL number A = At/ min(Ar, Az) we use is A = 0.5. 
Fourth-order Kreiss-Oliger dissipation [57] is added to the right-hand sides of the 
evolution equations, with a typical parameter value of e = 0.5. 

The elliptic equations are solved using a multigrid method [19]. The Full 
Approximation Storage (FAS) variant of the method enables us to solve nonlinear 
equations directly, i.e. we do not apply an outer Newton-Raphson iteration in order to 
obtain a sequence of linear problems. In the relaxation step of the multigrid algorithm, 
a nonlinear Gauss-Seidel relaxation (in red-black ordering) is directly applied to the full 
nonlinear equations. At each grid point, we solve simultaneously for the unknowns 7^^, 
V', a, and (3^. Only the Hamiltonian constraint ((23|) requires the solution of a (scalar) 
nonlinear equation, and this is done using Newton's method; a single iteration is found 
to be sufficient. All interior grid points are relaxed and afterwards the values at the 
ghost points are filled according to the boundary conditions. In order to transfer the 
numerical solution between the grids, we use biquadratic interpolation for prolongation 
and linear averaging for restriction. 

For the prolate wave evolved in section 15. 3|. the elliptic equations become highly 
anisotropic and the standard pointwise Gauss-Seidel relaxation employed in the 
multigrid method no longer converges. A common cure to this problem is line 
relaxation [28]. We solve for the unknowns at all grid points in a line z = const 
simultaneously. One Newton-Raphson step is applied to treat the nonlinearity and 
the resulting tridiagonal linear system is solved using the Thomas algorithm. Note 
that this method has the same computational complexity as the pointwise Gauss-Seidel 
relaxation. 

The wide range of length scales in the solutions we are interested in necessitates 
a position-dependent grid resolution. The classic adaptive mesh refinement (AMR) 
algorithm by Berger and Oliger [25j was designed for hyperbolic equations. Including 
elliptic equations in this approach is rather complicated. A solution with numerical 
relativity applications in mind was suggested by Pretorius and Choptuik [29j , and we 
shall use their algorithm here, with minor modifications due to the fact that our grids 
are cell centred rather than vertex centred. The key idea of the algorithm is that 
solution of the elliptic equations on coarse grids is deferred until all finer grids have 
reached the same time; meanwhile the elliptic unknowns are linearly extrapolated in 
time and only the evolution equations are solved. We have found that this approach 
works well as long as no grid boundaries are placed in the highly nonlinear region. 
In particular, adaptive generation of finer grids in the course of the evolution causes 
small but noticeable reflections that from our experience make the study of problems 
such as Brill wave critical collapse unfeasible. For this reason, the evolutions presented 
in this article use fixed mesh refinement (FMR), i.e. the grid hierarchy is defined at 
the beginning of the simulation and remains unchanged as time evolves. 

Finally we briefly discuss how an apparent horizon is found 'mat — const slice. 
The horizon is parametrized as a curve R = f{9) in spherical polar coordinates. 
Requiring the expansion of the outgoing null rays emanating from the horizon to 
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vanish yields a second order ordinary differential equation, which is solved using the 
shooting method. The boundary conditions are /'(O) = /'(7r/2) = 0, i.e., the horizon 
has no cusps on the axes. We follow an idea of Garfinkle and Duncan [Sj in order 
to monitor the approach to apparent horizon formation. For each point on the axis 
r ^ ■i=> 6 ^ 0, we find the angle 7 at which the curve starting from that point meets 
the z = 4^ 6 = 7r/2 axis. 



/' 

cot 7 = y 



(37) 



t/2 

We find the maximum 7max of this angle over all such curves. Obviously 7,nax = 7r/2 
for an apparent horizon, and the deviation from this value indicates how close we are 
to the formation of an apparent horizon. 

5. Numerical results 

As an application of our numerical implementation, we consider vacuum axisymmetric 
gravitational waves, so-called Brill waves The initial slice is taken to be time- 

symmetric so that W = r]"^ = initially. We consider the same initial data for the 
function S as in jTB] and [3], 

/ ^2 ^2 \ 

S = arcxp ] , (38) 

where a, cr^ and cr^ are constants. The initial lapse and shift are taken to be a = 1 
and 13^ = 0. The momentum constraints ([22|) are trivially satisfied initially and only 
the Hamiltonian constraint (l23l) needs to be solved. 



5.1. Convergence test 

In order to check convergence of the code, we first consider a wave with parameters 
a ~ ar ~ (Jz = 1. This will disperse rather than collapse to a black hole but is still 
well in the nonlinear regime. The ADM mass is Madm = 0.034. We take the domain 
size to be rmax = ^max = 40. The FMR hierarchy contains three grids (figure [1]). All 
the grids contain the origin, are successively refined by a factor of 2, and all have the 
same number of grid points iV^ — Nz. 

We run the simulation with three different resolutions, Nr — Nz ^ {64, 128, 256}. 
This enables us to carry out a three-grid convergence test: for each variable u we 
define a convergence factor 

^ _ II UAH - U2h \\ ,„„x 
II U2h - Uh I 

with the indices referring to the grid spacing. The norms are discrete norms taken 
over the subsets of all grids in the FMR hierarchy that do not overlap with finer grids. 
For a second-order accurate numerical method we expect Qu = 4. Figure [2] confirms 
that the code is approximately second-order convergent. (Occasional values Q„ > 4 
are not uncommon in similar numerical schemes [4ll29|.) 

As noted earlier, there are additional evolution equations for the variables 77+ , 
77_ and tp that are not actively evolved in our constrained evolution scheme. We use 
these to check the accuracy of the numerical implementation in the following way. We 
keep a set of auxiliary variables 77+ , fj- and '0 which are copied from their unhatted 
counterparts initially but evolved using the evolution equations (jA.3p - (jA.5p . During 
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Figure 1. FMR grid hierarchies used lor the Brill wave evolutions presented in 
this paper. Top left: err = cr^ = a = 1 f section 15.11 1: top right: (Jr = cr^ = 1, 
a = 8.5 (section l5.2ll . bottom: cr^ = 0.128, az = 1.6, a = 325 (section l5.3l l. 



the evolution, we form the differences between the two sets. Doing so for two different 
resolutions (grid spacings h and 2h) allows us to define another convergence factor for 
each u € {ry_|-, 7/_, -0} (referred to as residual convergence in the following), 



U2h - U2h 



(40) 



Uh - Uh I 

The results in figure [3] are again compatible with second-order convergence. 

We note that the residual convergence test just presented is more severe than the 
three-grid convergence test in the following sense. For the residuals of the unsolved 
evolution equations to converge as desired, not only must the numerical solution be 
second-order convergent but the constraint and evolution equations and their boundary 
conditions must be compatible. No exact boundary conditions are known at a finite 
distance from the source and compatibility of the boundary conditions we use is 
only achieved at infinity. We deliberately chose the domain size in this convergence 
test to be sufficiently large (~ 10'^Madm) so that the effect of the boundary on the 





Figure 2. Three-grid convergence factors I I39I I for a Brill wave with a 
(Tz = 1 computed from the three resolutions Nr = Nz = 64, 128, 256. 



convergence factors qu is small. 

As another consistency test, we compute a numerical approximation to the ADM 
mass. This is evaluated as a surface integral on a sphere close to the outer boundary, 
at spherical radius R = Rm , 

Mabm^-Rm Jau^ sine d0, (41) 
Jo 

where 

J A - 24,^A + ^rS,A, (42) 

ua is the unit normal in the spherical radial direction, and these expressions are valid 
in linearized theory. We evaluate Madm in (HJ) for two radii Ri\j £ {14,18} and 
extrapolate to infinity assuming that Madm(-Rm) — A^adm(oo) + const/i?j\/ . The 
result in figure [3] shows how numerical conservation of the ADM mass improves with 
increasing resolution. 



5.2. "Spherical" collapse 

Next we consider a wave with (7^ = CTz = 1 and a — 8.5. We refer to this as "spherical" 
because ar = ctz, although of course the actual wave is not spherically symmetric. 
Unlike the one in section lSTTl this wave is super-critical and will collapse to form a black 
hole. The ADM mass is -Madm ~ 2. We run the simulation for two different domain 
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Figure 3. Residual convergence factors l|40|l for a Brill wave with a = ar = cr^ = 
1 computed from two pairs of resolutions, Nr = = 64 : 128 (dashed) and 
128 : 256 (solid). The bottom right panel shows the numerically computed ADM 
mass for the three different resolutions, A'^^ = = 64 (dotted), 128 (dashed) 
and 256 (solid). 

sizes, rmax — Zmax G {10,20}. The FMR hierarchy is of a similar type as in section 
15.11 On the smaller domain there are three grids and for the larger domain we add on 
another coarse grid (figure [T]). We run the simulation with two different resolutions, 
Nr — Nz G {128,256}. In [6], the same initial data were evolved on a non-uniform 
grid with spacing Ar = Az = 1.92 x 10~^ close to the origin. This is comparable to 
our lower resolution grid hierarchy, which has grid spacing Ar = Az — 1.95 x 10~^ on 
the finest grid. 

Figure [4] shows the residual convergence factors defined in (|40|) . The general trend 
is that the convergence factors are close to 4 at late times but somewhat smaller at early 
times. Moving the outer boundary further out improves convergence considerably 
at early times, as can be seen particularly for the variables ri±. This demonstrates 
the effect of the outer boundary where imperfect boundary conditions are imposed 
at a finite distance. Because of the elliptic equations involved in our evolution 
scheme, inaccuracies in the outer boundary conditions influence the entire domain 
instantaneously, not only after the outgoing radiation interacts with the boundary as 
is the case in a hyperbolic scheme. Moving the outer boundary much further out by 
adding more coarse grids is not feasible for this evolution because of the computational 
cost involved in the current single-processor implementation of the code. 

In particular, the value of the conformal factor ip far out appears to be very 




sensitive to the dynamics close to the origin. This is not the case for ip, which is 
evolved from the initial data by a hyperbolic PDE. As a result, the difference ip — ip has 
a large, spatially nearly constant contribution that is nearly resolution independent, 
thus causing the convergence factor to degrade. 

At times later than those shown here, the convergence factors ultimately decrease 
because large gradients develop due to the grid-stretching property of maximal slicing. 
However, here we are only interested in the part of the evolution until just after the 
formation of the apparent horizon. 

Figure |4] also indicates that both increasing the resolution and the boundary 
radius improves the numerical conservation of the ADM mass. For the larger domain 
at the higher resolution, the initial oscillations are at the 5% level and after < « 1 the 
mass remains constant to within 1%. 

Next we evaluate the lapse function a in the origin r ^ z — 0. As a consequence of 
the singularity avoidance property of maximal slicing, the lapse is expected to collapse 
towards zero as a strong-gravity region of spacetime is approached. Our result in figure 
([5|) is in good agreement with ^ and appears to be insensitive to the resolution and 
boundary location. 

We also plot the Riemann invariant / = R"-^'^'^ Rabcd / in the origin. The decay 
of this quantity after i w 1 agrees roughly with [B], although we find a somewhat 
different behaviour at earlier times (rather than increasing right from the beginning. 




Figure 6. Apparent horizon finder angle 7max (cf. I I37I I) and mass Mah for a 
Brill wave with ar = cr^ = 1 and a = 8.5. Results for two different domain sizes 
''max = ^max and for two different resolutions A''^ = Nz = 128 (dashed) and 256 
(solid) are shown. 



/ first decreases for a short time). However there is a rather strong dependence on 
resolution and outer boundary location in this case, which indicates that the results 
for / should be interpreted with care. 

Finally we search for an apparent horizon. The evolution of the angle 7max 
(cf. ([571) ') is shown in figure [SI It agrees reasonably well with although we find that 
the horizon forms slightly earlier at t = 3.6 ± 0.2 rather than at t = 3.9. Also shown 
in figure[niis the mass of the horizon, computed from its area Aah = IGttM^jj. When 
it first forms, the horizon has mass Mah = 1-85 ± 0.05. The numerically computed 
ADM mass at this time is Madm = 2.04 ± 0.02 so that Mah/Madm = 0.91 ± 0.03, as 
compared with Mah/Madm = 0.82 reported in [6]. After its formation the horizon 
expands slightly (its mass increases by about 3%) and appears to ultimately settle 
down. The results stated here correspond to the run on the larger domain at the 
higher resolution and the errors are estimated by comparison with the other runs. 
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Figure 7. Lapse function a and Riemann invariant / in the origin for a Brill 
wave with cr,. = 0.128, (Tz =1.6 and a = 325. 



5.3. Highly prolate collapse 

We now turn to a highly prolate Brill wave with ar = 0.128, (Tz = 1.6, and a = 325, 
which again has Madm ~ 2. This is one of the initial data sets considered in (TB] and 
it was also evolved (until t « 1.5) in [6]. 

Our spatial domain has size rmax — ^max = 20. The resolution on the base grid 
is taken to be Nr = 256 and — 64. There are 6 grids in the FMR hierarchy. Again 
all the grids contain the origin and the grid spacing is successively halved in both 
dimensions. The number of grid points Nr in the radial direction is the same on all 
grids but is successively multiplied by a factor of (approximately) 1.34. In this 
way the finer grids are better adapted to the prolate shape of the initial data. The 
grid hierarchy is shown in figure [T] The spacing on the finest grid is Ar = 2.44 x 10~^ 
and Az ~ 9.77 x 10^'^. By comparison, the grid used in [B] has Ar = 9.70 x 10~^ and 
Az = 3.74 X 10^^ close to the origin, roughly four times coarser. 

Figure [7] shows the evolution of the lapse function a and Riemann invariant / 
in the origin. These agree well with [6], except for the t < 0.5 part of / (but note 
the strong dependence of this quantity on resolution and outer boundary location 
apparent from figure [5|). 

The approach to apparent horizon formation is shown in figure [H We are able to 
evolve the wave for much longer than the authors of [B] and we confirm their conjecture 
that an apparent horizon will indeed form. It first appears at t = 3.9 and its shape is 
remarkably close to a sphere in our coordinates. At its formation the apparent horizon 
has mass A/ah = 1-990 and at this time the ADM mass has settled down to a value of 
-^ADM = 2.065 so that Mah/A/adm = 0.96. This is in accordance with the Penrose 
inequality [3Dj, which conjectures that AfAn/AfADM ^ 1- 

6. Conclusions 

We considered constrained evolution schemes for the Einstein equations in 
axisymmetric vacuum spacetimes. One of the motivations for this work was to 
try and understand why the numerical elliptic solvers in some of these schemes, 
e.g. [4], failed to converge in certain situations. We found that this was related to 
the elliptic equations becoming indefinite. Apart from the implications for numerical 
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Figure 8. Apparent horizon formation for a Brill wave with (Jr = 0.128, cr^ = 1.6 
and a = 325. Top left: horizon finder angle 7max (cf. I|37|l ). Top right: ADM mass 
(solid) and apparent horizon mass (dashed). Bottom: coordinate location of the 
apparent horizon when it first forms at t = 3.9 (solid), and at t = 4.5 (dashed). 

convergence, we also pointed out that such equations might admit nonunique solutions. 
In section l373l we presented a new scheme that does not suffer from this problem. Its 
main features are a suitable rescaling of the extrinsic curvature with the conformal 
factor, and separate solution of the momentum constraints and isotropic spatial gauge 
conditions. Thus the scheme involves the solution of six elliptic equations rather than 
four as in [5 . Given that multigrid methods [12] can be used to solve these equations 
at linear complexity, this does not imply a severe increase in computational cost. 

Our numerical implementation uses second-order accurate finite differences and 
combines mesh refinement with a multigrid elliptic solver, based on the algorithm of 
[29j . We work in cylindrical polar coordinates. Unlike in we do not compactify the 
spatial domain but impose boundary conditions at a finite distance from the origin. 

As an application of the code, we evolved Brill waves in section [5l We carried 
out a careful convergence test in section 15.11 and demonstrated that the code is 
approximately second-order convergent. For a stronger super-critical wave (section 
15. 2|) . convergence of the residuals of the unsolved evolution equations was somewhat 
slower at earlier times. Varying the domain size indicated that this is mainly caused 
by inaccuracies in the outer boundary conditions we use. These errors appear to have 
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little effect on tfie formation of the apparent horizon. 

In section 15.31 we evolved a highly prolate Brill wave. Such initial data were 
conjectured in [16] to form a naked singularity rather than a black hole, which 
would violate weak cosmic censorship. However an apparent horizon does form in 
our evolution. We thus improve on the results of the authors of [6], who could not 
evolve the wave for a sufficiently long time to see an apparent horizon form although 
they conjectured that this would happen eventually. 

There are many directions in which this work could be extended. For simplicity 
we only considered vacuum spacetimes with a hypersurface-orthogonal Killing vector, 
i.e., vanishing twist. The addition of matter and twist should be straightforward. Care 
must be taken that the additional variables are rescaled with suitable powers of the 
conformal factor so that the Hamiltonian constraint (j23|) remains definite. An elegant 
framework capable of including twist is provided by the (2 + 1) + 1 formalism [3T]l32j. 
From a mathematical point of view, it would be interesting to prove that the Cauchy 
problem or even the initial-boundary value problem is well posed for the present (or 
a similar) formulation of the axisymmetric Einstein equations. These questions were 
studied for similar hyperbolic-elliptic systems in [33l [SU [35] . 

It is a disadvantage of constrained evolution schemes that inaccuracies in the 
outer boundary conditions influence the entire domain instantaneously. More work is 
needed on improved boundary conditions in the context of mixed hyperbolic-elliptic 
formulations of the Einstein equations. An alternative to an outer boundary at a 
finite distance would be the compactification towards spatial infinity used in [6]. 
However, outgoing waves ultimately fail to be resolved on such a compactified grid, 
which because of the elliptic equations involved has again adverse effects on the entire 
solution. This problem is avoided if hyperboloidal slices are used, which can be 
compactified towards future null infinity (see [36] for a related review article). In this 
case the constraint and evolution equations become formally singular at the boundary, 
which needs to be addressed in a numerical implementation. 

On the computational side, the accuracy of the code could be improved by using 
fourth (or higher) order finite differences. Computational speed could be gained by 
parallelizing the code and running it on multiple processors. 

It would be interesting to evolve even more prolate Brill waves than the 
one considered here. However the elliptic equations then become more and more 
anisotropic and the relaxation method employed in the multigrid method must be 
modified to ensure convergence. For the wave considered in this paper, line relaxation 
was found to accomplish this but we have not been able to achieve convergence for 
even more prolate configurations. More sophisticated modifications such as operator- 
based prolongation and restriction [28j are likely to be required. In any case, in order 
to evolve some of the extremely prolate initial data sets considered in [16j where 
ar/cTz ~ 10"^ in ([38| . a radically different approach will probably be needed. 

Another interesting application of our code would be Brill wave critical collapse. 
However, preliminary results indicate that we are currently unable to evolve waves 
close to the critical point for a sufficiently long time because reflections from the 
interior AMR grid boundaries become increasingly pronounced as more and more 
finer grids are added close to the origin. The mesh refinement algorithm of that 
we adopt appeared to be sufficiently robust in the scalar field evolutions of but 
we suspect that the situation is quite different in vacuum collapse. An improved AMR 
algorithm for mixed hyperbolic-elliptic systems of PDEs will probably be required. 
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Appendix A. Evolution equations 

Here we list the evolution equations of the formulation of the axisymmetric Einstein 
equations presented in section [3?3l The variables S and W are actively evolved, 

S^t = P^'S^r + f3'S,^ - aip-^W + r-^p^'S + (r" V^),r , (A.l) 

-\aip~^r-^[ff+ + 2f^^rW - A{rWf] + Ar'^P'^W. (A.2) 

The variables -0 and t]^ are not evolved but solved for using the constraint equations. 
However the Einstein equations also imply the following evolution equations that may 
be used to check the accuracy of a numerical scheme. 

■tP ^ = p'-^p^^ + + ^r-i/?'~V - - 2rW), (A.3) 

77+,t = 2aip^e-^''^[~a-^a^rz ~ 2?A~V,rz + RrA^ + Rz{Ar + r^^) + QPrPz 

+ 2Pr{A, +R,)+ 2P,{Ar + Rr)] 

+P''{ri\rz + V'.rr + ir-^V+) + P" W ,zz + V\rz) " V+P- 

+ |(/3'',^ - 2P\r)v- + lrWP+ + \aij-^7]+{r^- + ArW), (A.4) 
7y_ t = aV^e"^'~^[-a"^a^.rr + a^^a,^^ - 2V'~V,rr + 2^'" + 2Rr\Ar + r"^) 
-2A,^R, + 2Pr{2Ar + iPr + 2Rr) ~ 2P^{2A^ + 3P^ + 2R^)] 

-f/3_(27y_ - rW) + {P\r - P'',z)v+ + ^ai/'"^r7-(77- + ArW). (A.5) 

Appendix B. Finite-difference operators and boundary conditions 

Centred second-order accurate finite-difference operators are used at all interior grid 
points. We only give the expressions for derivatives in the r direction; corresponding 
expressions hold in the z direction. The symbol « denotes equality up to 0(Ar)^. 

u,r|i ~ {ui+1 ~ Ui-i)/(2Ar), 

{r~^u)^r\i ~ [ui+i/ri+i - Ui_i/ri_i)/(2Ar), 

U.rrli « - 2Ui + Ui_i)/(Ar)^, 

[r~'^u^r)Ai ~ - Ui)/ri+i - {ui - Ui_2)/?'i-i]/[4(Ar)^]. (B.l) 

At r = (z = 0), a Dirichlet condition is imposed if the variable is an odd function 
of r (z) and a Neumann condition if it is an even function of r (z). The parities of 
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the various variables are as follows, 

odd in r: , tj'' , S , W 
even in r : a, -0, , ry^ 
odd in z : (3^',r]^' 

evenin z : a,ip,(3\T]'',S,W (B.2) 

The value of a variable u at the ghost point is set to be uo = —ui if u obeys a Dirichlet 
condition and uq = ui if it obeys a Neumann conditiorlj. This discretization at the 
boundary is second-order accurate. 

Dirichlet conditions at the outer boundary are implemented in a similar way. In 
order to impose the falloff condition ((35)) . we note that Ru is a linear function of R 
and so we linearly extrapolate Ru in the radial direction from the interior grid points 
to the ghost points in order to find the values of u there. The Sommerfeld condition 
(|36p is rewritten in the form 

and discretized at the outer ghost points on the base grid of the AMR hierarchy in 
order to integrate the values of u there forward in time. Here backward differencing 
is used in the direction normal to the boundary, 

U^rV ~ {^iu^ - AUr-i + Uj_2)/(2Ar) (B.4) 

and similarly in the z direction. 
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